function [averageBetas maskedBetas averageBetasMatrix stdBetasMatrix] = extractROImask(subjID, subjectsDir)
%
%   _______________________________________________
%   by Marcelo G Mattar (12/20/2012)
%   mattar@sas.upenn.edu


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CHECK INPUTS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin < 2
    subjectsDir = '/Users/marcelomattar/Data/ITIanalysis/Subjects/';
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BASIC DEFINITIONS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% Define some basic directories and files
ROIsdir = [subjectsDir subjID '/ROIs/'];
dataDir = [subjectsDir subjID '/BinarySeq/'];
featDirs = dir([dataDir '*.feat']);
indexOf_1back_adaptation_within_feat = 3;

numFeatDirs = length(featDirs);
averageBetas = cell(numFeatDirs+1,5);
averageBetas{1,1} = 'ITI';
averageBetas{1,2} = 'Covariates';
averageBetas{1,3} = 'Temp-Deriv?';
averageBetas{1,4} = 'lFFA';
averageBetas{1,5} = 'rFFA';
averageBetas{1,6} = 'Both FFAs';
averageBetas{1,7} = 'Standard Deviation';

% Pre-allocate some important variables
maskedBetas = cell(numFeatDirs,3);
averageBetasMatrix = zeros(numFeatDirs,3);
stdBetasMatrix = zeros(numFeatDirs,3);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD FFA MASK
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load the masks
lFFAmask = loadniigz([ROIsdir 'lFFA_FvsS.nii.gz']);
rFFAmask = loadniigz([ROIsdir 'rFFA_FvsS.nii.gz']);

FFAmask = lFFAmask + rFFAmask;

% Make each mask into a logical matrix
FFAmask = logical(FFAmask);
lFFAmask = logical(lFFAmask);
rFFAmask = logical(rFFAmask);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% LOAD ADAPT DATA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Extract Adaptation copes, varcopes and tstats for each of the Face runs
for thisDir = 1:numFeatDirs
    feat_cope = loadniigz([dataDir featDirs(thisDir).name '/stats/cope' num2str(indexOf_1back_adaptation_within_feat) '.nii.gz']);
    maskedBetas{thisDir,1} = feat_cope(lFFAmask(:));
    maskedBetas{thisDir,2} = feat_cope(rFFAmask(:));
    maskedBetas{thisDir,3} = feat_cope(FFAmask(:));
    
    averageBetasMatrix(thisDir,1) = mean(maskedBetas{thisDir,1});
    averageBetasMatrix(thisDir,2) = mean(maskedBetas{thisDir,2});
    averageBetasMatrix(thisDir,3) = mean(maskedBetas{thisDir,3});
    
    stdBetasMatrix(thisDir,1) = std(maskedBetas{thisDir,1});
    stdBetasMatrix(thisDir,2) = std(maskedBetas{thisDir,2});
    stdBetasMatrix(thisDir,3) = std(maskedBetas{thisDir,3});
    
    averageBetas{thisDir+1,4} = mean(maskedBetas{thisDir,1});
    averageBetas{thisDir+1,5} = mean(maskedBetas{thisDir,2});
    averageBetas{thisDir+1,6} = mean(maskedBetas{thisDir,3});
    averageBetas{thisDir+1,7} = std(maskedBetas{thisDir,3});
    
    [prefix, remain] = strtok(featDirs(thisDir).name, '_');
    [title, remain] = strtok(remain, '_');
    ITI = title(10:end);
    type = strtok(remain, '.');
    type = type(2:end);
    covs = strtok(type,'b');
    averageBetas{thisDir+1,1} = ITI;
    averageBetas{thisDir+1,2} = covs;
    averageBetas{thisDir+1,3} = (prefix(end)=='d');
end

display('end');
